﻿using LatoolNet;

namespace NumericalSolution {
  class ConvectionDiffusion1DQuick : ConvectionDiffusion1D {
    public override Matrix Solve() {

      int subDiagonal = 2;
      int superDiagonal = 1;
      Matrix mat = new Matrix(p_numberOfControlVolumes,
                              p_numberOfControlVolumes,
                              subDiagonal,
                              superDiagonal);
      Matrix vec = new Matrix(p_numberOfControlVolumes, 1);

      double phi_A = 1.0;
      double phi_B = 0.0;

      double phi_WW;
      double phi_W;
      double phi_P;
      double phi_E;
      double c;

      phi_P = (7.0 / 8.0) * p_rho * p_u + 4.0 * p_gamma / p_deltax;
      phi_E = (3.0 / 8.0) * p_rho * p_u - (4.0 / 3.0) * p_gamma / p_deltax;
      c = (10.0 / 8.0) * p_rho * p_u * phi_A 
        + (8.0 / 3.0) * p_gamma * phi_A / p_deltax;

      mat[0, 0] = phi_P;
      mat[0, 1] = phi_E;
      vec[0, 0] = c;

      phi_W = -p_rho * p_u - p_gamma / p_deltax;
      phi_P = (3.0 / 8.0) * p_rho * p_u + 2.0 * p_gamma / p_deltax;
      phi_E = (3.0 / 8.0) * p_rho * p_u - p_gamma / p_deltax;
      c = -(2.0 / 8.0) * p_rho * p_u * phi_A;

      mat[1, 0] = phi_W;
      mat[1, 1] = phi_P;
      mat[1, 2] = phi_E;
      vec[1, 0] = c;

      for (int i = 2; i < p_numberOfControlVolumes - 1; i++) {

        phi_WW = (1.0 / 8.0) * p_rho * p_u;
        phi_W = -(7.0 / 8.0) * p_rho * p_u - p_gamma / p_deltax;
        phi_P = (3.0 / 8.0) * p_rho * p_u + 2.0 * p_gamma / p_deltax;
        phi_E = (3.0 / 8.0) * p_rho * p_u - p_gamma / p_deltax;
        c = 0.0;

        mat[i, i - 2] = phi_WW;
        mat[i, i - 1] = phi_W;
        mat[i, i] = phi_P;
        mat[i, i + 1] = phi_E;
        vec[i, 0] = 0;
      }

      phi_WW = (1.0 / 8.0) * p_rho * p_u;
      phi_W = -(6.0 / 8.0) * p_rho * p_u - (4.0 / 3.0) * p_gamma / p_deltax;
      phi_P = -(3.0 / 8.0) * p_rho * p_u + 4.0 * p_gamma / p_deltax;
      c = (8.0 / 3.0) * p_gamma * phi_B / p_deltax - p_rho * p_u * phi_B;

      mat[p_numberOfControlVolumes - 1, p_numberOfControlVolumes - 3] 
        = phi_WW;
      mat[p_numberOfControlVolumes - 1, p_numberOfControlVolumes - 2] 
        = phi_W;
      mat[p_numberOfControlVolumes - 1, p_numberOfControlVolumes - 1] 
        = phi_P;
      vec[p_numberOfControlVolumes - 1, 0] = c;

      //LUFactorization.Solve(mat, vec);
      LinearEquation.Solve(mat, vec);

      return vec;
    }

  }
}
